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ABSTRACT 

We use two-dimensional time-dependent hydrodynamical simulations to follow the growth of the 
Kelvin-Helmholtz (K-H) instability in cooling jets into the nonlinear regime. We focus primarily on 
asymmetric modes that give rise to transverse displacements of the jet beam. A variety of Mach numbers 
and two different cooling curves are studied. The growth rates of waves in the linear regime measured 
from the numerical simulations are in excellent agreement with the predictions of the linear stability 
analysis presented in the first paper in this series. In the nonlinear regime, the simulations show that 
asymmetric modes of the K-H instability can affect the structure and evolution of cooling jets in a 
number of ways. We find that jets in which the growth rate of the sinusoidal surface wave has a 
maximum at a so-called resonant frequency can be dominated by large-amplitude sinusoidal oscillations 
near this frequency. Eventually, growth of this wave can disrupt the jet. On the other hand, nonlinear 
body waves tend to produce low-amplitude wiggles in the shape of the jet but can result in strong 
shocks in the jet beam. In cooling jets, these shocks can produce dense knots and filaments of cooling 
gas within the jet. Ripples in the surface of the jet beam caused by both surface and body waves generate 
oblique shock “spurs” driven into the ambient gas. Our simulations show these shock “spurs” can 
accelerate ambient gas at large distances from the jet beam to low velocities, which represents a new 
mechanism by which low-velocity bipolar outflows may be driven by high-velocity jets. Rapid entrain- 
ment and acceleration of ambient gas may also occur if the jet is disrupted. 

For parameters typical of protostellar jets, the frequency at which K-H growth is a maximum (or 
highest frequency to which the entire jet can respond dynamically) will be associated with perturbations 
with a period of ~200 yr. Higher frequency (shorter period) perturbations excite waves associated with 
body modes that produce internal shocks and only small-amplitude wiggles within the jet. The fact that 
most observed systems show no evidence for large-amplitude sinusoidal oscillation leading to disruption 
is indicative that the perturbation frequencies are generally large, consistent with the suggestion that pro- 
tostellar jets arise from the inner regions (r < 1 AU) of accretion disks. 

Subject headings: galaxies : jets — hydrodynamics — instabilities — ISM: jets and outflows — MHD 


1. INTRODUCTION 

Some of the best studied examples of astrophysical jets 
are observed at optical wavelengths. Examples include both 
protostellar (Herbig-Haro) jets (Edwards, Ray, & Mundt 
1993; Ray 1996) and jets from some types of active galactic 
nuclei, e.g., Seyfert galaxies. Optical jets make excellent 
laboratories for studying the hydrodynamics of collimated 
outflows for several reasons: (1) they are amenable to high- 
resolution imaging both by ground-based adaptive optics 
systems and the space-based Hubble Space Telescope 
( HST ), and (2) in most cases, the optical radiation arises 
from forbidden emission lines, which allows spectroscopic 
studies to be used to determine both the kinematics and the 
physical properties (e.g., density, temperature, etc.) of the 
line-emitting gas. 

Perhaps the most striking images of optical jets are pro- 
vided by recent HST observations of several protostellar 
systems. With the high spatial resolution of HST, structures 
within the jet are resolved, which allows direct study of the 
internal dynamics of the jet beam in addition to studies of 
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the effect of the jet on its surrounding ambient gas. In three 
systems, HH 30 (Stapelfeldt et al. 1997; Ray et al. 1997), HH 
111 (Reipurth et al. 1997), and HH 34, the jet appears as a 
continuous but knotty band of emission near the source. In 
each system, however, the jet is neither perfectly straight 
(e.g., shows low-amplitude long-wavelength sinusoidal 
oscillations) nor perfectly symmetric (e.g., contains knots 
and filaments of emission that are not distributed symmetri- 
cally with respect to the center of the jet). HST images of the 
HH 47 jet (Heathcote et al. 1996) reveal an even more 
complex flow in this system. The jet is very filamentary and 
appears to undergo large-amplitude bends before termina- 
ting in a classic Mach disk and bow shock structure at the 
head of the jet. Associated with the kinks and bends in the 
HH 47 jet are Ha-emitting shock spurs driven into the 
ambient gas. 

While it is obvious that images of more distant extra- 
galactic jets will not reveal as much detail as those of nearby 
protostellar jets, recent observations of Seyfert jets have 
discovered complex and puzzling asymmetric structures. 
Perhaps the best example is the helical jet in NGC 4258 
(Cecil, Wilson, & Tully 1992) in which three tightly wrapped 
strands appear to spiral around one another for two full 
cycles before dispersing. Kinematics of the jet suggest that 
the jet material moves along the strands, so that the helical 
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pattern is not a fixed structure carried outward by the flow. 
These observations led Cecil et al. to speculate that fluid 
dynamical instabilities in the shear layer between the jet and 
ambient gas might be responsible for forming the strands. 

The observations suggest several possibilities for the for- 
mation of complex structure in jets: 

1. A possibility suggested by both kinematic and mor- 
phological data of protostellar jets is that they are time 
variable in velocity, perhaps owing to variability in the 
underlying driving mechanism. Fluctuations in the jet 
velocity that exceed the internal sound speed result in the 
formation of internal shocks (Raga & Kofman 1992), and 
hydrodynamic simulations of the propagation of pulsed jets 
show that these shocks develop quickly near the jet nozzle 
and produce a layer of postshock cooling gas in the jet that 
can account for many of the observed properties of the 
observed emission knots (Hartigan & Raymond 1993; 
Stone & Norman 1993b; Gouviea dal Pino & Benz 1994; 
Biro & Raga 1994). 

2. Inhomogeneities in the external medium can affect jet 
propagation and produce complex structure. For example, 
Stone & Norman (1994) investigated the effect of a trans- 
verse density gradient on the propagation of a cooling jet, 
and, more recently, Raga & Canto (1996) have considered 
the effect of a jet/cloud impact with possible application to 
the HH 110 system (Reipurth, Raga, & Heathcote 1996). 
This work was preceded by even earlier work that examined 
the interaction of extragalactic jets with clouds in the 
interstellar/intergalactic medium (Wilson 1988; Norman & 
Balsara 1993 ; Clarke 1993). 

3. It has been suggested that precession of the accretion 
disk can result in large reorientations in the direction of 
propagation of a jet, which produces either discrete knots 
that follow diverging paths (Raga & Biro 1993) or a more 
continuous oscillatory pattern (Raga, Canto, & Biro 1993; 
Biro, Raga, & Canto 1995). 

4. Finally, as first suggested by Biihrke, Mundt, & Ray 
(1988; hereafter BMR), the growth of Kelvin-Helmholtz 
(K-H) instabilities in the jet can result in the formation of 
structure (see also Bodo et al. 1994; Massaglia et al. 1992). 
BMR suggested that waves associated with the pinch mode 
of the K-H instability may result in knots in the jet, 
although it now appears that observations of the spacing 
and proper motion (see, e.g., Eisloffel & Mundt 1994) of the 
observed knots favors the internal working surface interpre- 
tation. Nevertheless, extensive studies of the K-H instability 
relevant to extragalactic jets (see Birkinshaw 1991 for a 
review) reveals that supermagnetosonic jets are always 
unstable. Thus, asymmetric modes of the K-H instability 
may still play a significant role in the development of struc- 
ture. 

An important property of optical jets is that in most 
systems (e.g., both protostellar and optical Seyfert jets), the 
optical emission results in a loss of internal energy from the 
gas, i.e., the jet is radiatively cooled. If the cooling time of 
the gas is comparable to or less than the dynamical flow 
time, the loss of internal energy via cooling can substan- 
tially alter the dynamics of cooling jets in comparison to 
adiabatic jets. Thus, although the stability properties of 
“adiabatic” extragalactic jets are relatively well studied, 
such is not the case for cooling jets. Both Hunter & Whit- 
aker (1989) and Massaglia et al. (1992) showed that radi- 
ative cooling can significantly modify the K-H instability. 


In a companion paper (Hardee & Stone 1997; hereafter 
Paper I), we presented a comprehensive linear stability 
analysis of “cooling” jets. In Paper I, the growth of K-H 
symmetric pinch and asymmetric sinusoidal modes were 
calculated in the linear regime for several Mach numbers 
over a broad range of perturbation frequencies. Not sur- 
prisingly, we found that the growth rates and wavelengths 
of the K-H modes in the cooling jet differed substantially 
from those of the adiabatic jet. Depending on the slope of 
the cooling curve near the equilibrium temperature, modes 
could either be damped or amplified relative to the adia- 
batic limit. Fluid displacement surfaces associated with 
various surface and body waves associated with the pinch 
and sinusoidal modes were also calculated in order to 
predict structures in the nonlinear regime. This analysis 
indicated that low-frequency perturbations should be the 
most effective at disrupting the jet. In this paper, we extend 
the analysis of Paper I by following the growth of K-H 
unstable waves into the nonlinear regime using time- 
dependent hydrodynamical simulations. Our goal is to 
study both the saturation and structure of K-H modes in 
the nonlinear regime and to provide direct comparison to 
observations. 

There is an important difference between our present 
numerical studies of the K-H instability in jets and numeri- 
cal studies of “ precessing ” jets, even though the tools and 
techniques of both studies are otherwise very similar. First, 
it is important to emphasize that the K-H instability studied 
here is a fluid dynamical phenomenon that can be captured 
only with a full hydrodynamical treatment of the evolution 
of a jet. Despite the fact that the jets studied here are highly 
overdense with respect to their surroundings, their 
dynamics is not well approximated by ballistic particles (an 
approximation that is valid only for jets that are more than 
100 times denser than their surroundings). Moreover, in 
models of precessing jets, the direction of propagation of the 
jet is assumed to undergo large-amplitude variations, which 
result in transverse velocities at the jet nozzle that exceed 
the sound speed . Thus, the body of the jet cannot respond 
dynamically, and transverse shocks whose properties are 
determined by the assumed amplitude and frequency of the 
driving motion result. In the work presented here, we 
provide only small-amplitude, linear perturbations to the 
transverse motion of the jet. Typically, the transverse veloci- 
ties induced by such perturbations are only a few percent of 
the sound speed; thus, the perturbations will result in 
linear-amplitude transverse sound waves in the jet. These 
linear waves are then amplified by the shear at the jet- 
ambient medium interface, ultimately producing large- 
amplitude, nonlinear phenomena occurring far downstream 
of the nozzle, and whose properties reflect the dynamical 
properties of the resonant cavity that is formed by the jet 
beam. The difference between the hydrodynamical evolu- 
tion of nonlinear and linear amplitude perturbations 
applied at the jet nozzle is both substantial and essential. 
Moreover, it seems certain that astrophysical jets will be 
subject to linear amplitude perturbations owing to a variety 
of effects (nonlinear amplitude perturbations are less likely); 
thus, K-H modes should be relevant to all jets to some 
degree. Of course, in a real system, K-H instabilities may 
occur in addition to jet velocity variations (pulsing) and 
perhaps even large-scale readjustments to the jet direction. 
Thus, detailed models of specific observations may require 
the addition of more than one of these effects. 
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Finally, we note that for the reasons discussed in Paper I, 
we confine the analysis in this paper to two-dimensional 
slab jets. By adopting a slab geometry, we are able to study 
the asymmetric sinusoidal mode of the jet, as well as the 
symmetric pinch mode. While the pinch mode has been 
studied as a possible source of emission knots in protostel- 
lar jets (BMR; Massaglia et al. 1992), it is the sinusoidal 
mode that can result in wiggles, bends, and jet disruption. 
Moreover, the linear analysis shows that the sinusoidal 
mode on the slab jet is a direct analog of the helical mode 
on a three-dimensional jet of circular cross section. By con- 
fining our studies to two dimensions, we are able to perform 
a numerical parameter survey more efficiently. We will 
study the linear stability properties and nonlinear evolution 
of a three-dimensional cooling jet in a future communica- 
tion. 

The organization of this paper is as follows. In the follow- 
ing section, we describe our numerical methods for evolving 
the equations of hydrodynamics including optically thin 
radiative cooling. In § 3, we compare our numerical results 
with the analytic analyses of Paper I in the linear regime 
and describe the evolution of K-H unstable jets in the non- 
linear regime for a variety of parameters. In § 4, we compare 
our results to observations of protostellar jets, and in § 5, we 
conclude. 

2. NUMERICAL METHODS 
2. 1 . Hydrodynamical Algorithms 

The dynamical evolution of a supersonic jet is given by 
solutions to the mass, momentum, and total energy conser- 
vation equations : 

+ V • (/>») = 0 , (1) 

dpv 

+ V • (P + pvv) = 0 , (2) 

dpE 

- J ^ r + \-{Pv + pvE) = -n 2 A + nH , (3) 

where p is the mass density, v is the velocity vector, P is the 
pressure tensor with diagonal components p, and E = 
(y — 1 )~ 1 pjp + 0.5(t; * v) is the total specific energy (here, we 
use y = 5/3). The two terms on the right-hand side of 
equation (3) account for energy losses due to optically thin 
radiation and heating of the gas, respectively. In these 
terms, A is the per particle cooling rate, H is the heating 
rate, and n = p/m is the number density of particles of mean 
mass m. We adopt m = 1.4m H , where m n is the mass of the 
hydrogen atom. 

In order to generate solutions to equations (1M3), we use 
an implementation of the piecewise parabolic method 
(PPM) to evolve the left-hand side of equations (l)-(3) 
numerically. Following standard practice, we operator-split 
the solution of the right-hand side of equation (3) into a 
separate step. We use a fully implicit backward Euler 
scheme combined with Newton-Raphson iterations to 
advance the nonlinear heating and cooling terms. This 
ensures that our numerical scheme is stable even when the 
cooling time is much less than the dynamical time. The 
exact form adopted for A and H is described below. 

For highly supersonic flows, we have found the accuracy 
of the numerical evaluation of the heating and cooling rates 


is improved by solving an internal energy equation along 
with the other conservation laws. In hypersonic flows, the 
total energy E is dominated by the kinetic energy E k . The 
small difference between two large numbers (E — E k ) leads 
to round-off errors in the calculation of the temperature, 
which is proportional to E — E k . In most cases, e.g., studies 
of propagating jets, this error does not affect the dynamics 
since (1) it is extremely small (of order the machine word 
length) and (2) it occurs only where the pressure is negligi- 
ble. However, in this work, we wish to study the dynamical 
evolution of unstable jets that are initially in a delicate pres- 
sure equilibrium with a low-density ambient medium and in 
which the net cooling rate (which is a strongly nonlinear 
function of temperature) is n 2 A — nH = 0. In this case, 
round-off errors can introduce perturbations that are 
amplified by the net cooling rate and that can eventually 
overwhelm those introduced at the jet nozzle. Thus, we 
improve the consistency of the temperature evolution by 
solving an internal energy equation 

de « P „ 

— + v • \e V * v , (4) 

ot p 

using fluxes evaluated numerically from the PPM algorithm 
(Bryan et al. 1995), where e is the internal energy, and then 
use e to compute the temperature in the evaluation of the 
cooling terms. The total energy E is still used to evolve the 
hydrodynamics. Only with this extension were we able to 
hold the initial equilibrium state for many dynamical times 
if it was not perturbed. 

In order to follow the mixing of jet and ambient material, 
we evolve a Lagrangian tracer variable that we term the 
“ color ” C of the fluid. The evolution of C is described by 
the additional conservation law 

dpC 

-^ + V-pO> = 0, (5) 

which is updated along with equations (1)— (3) using PPM. 
Initially, the jet material is labeled by C = 1, while the 
ambient gas is labeled C = 0. As the jet evolves, the ambient 
and jet gases mix, which results in regions with 0 < C < 1. 

2.2. Cooling Rates 

An important aspect of the numerical simulation of the 
dynamical evolution of cooling jets is an accurate treatment 
of the microphysical heating and cooling rates. In a series of 
papers that studied the propagation of protostellar jets. 
Stone & Norman (1993a, 1993b, 1994) adopted a nonequi- 
librium ionization formalism to improve the accuracy of the 
cooling terms. This method required solving an ionization 
and recombination rate equation for neutral hydrogen in 
addition to the conservation laws equations (l)-(3) in order 
to compute the ionization fraction of the gas. Cooling due 
to both the ionized and neutral components of the medium 
could then be calculated self-consistently along with the 
dynamics. In Stone & Norman (1993a), it was shown 
through the use of a test problem based on the overs t ability 
of radiative shocks that this nonequilibrium treatment of 
the cooling rate led to significantly improved estimates for 
the total cooling rates in comparison to cooling based on 
the assumption of complete ionization. However, incorpor- 
ating a time-dependent ionization fraction into the linear 
stability analysis of a cooling jet is intractable. Thus, in 
Paper I, the cooling rate was assumed to be given by equi- 



No. 1, 1997 


RADIATIVELY COOLING JETS. II. 


139 


librium cooling curves in order to study the stability 
properties of the cooling jet. Two separate cooling curves 
were adopted in Paper I: one appropriate to interstellar gas, 
and one to photoionized gas with reduced metallicity. Since 
in Paper I we consider only linear amplitude perturbations 
in the cooling rate away from the equilibrium state, each of 
these cooling functions was represented by a power law in 
the neighborhood of the equilibrium temperature. 

In this paper, we adopt the same assumptions as used in 
Paper I in order to allow direct comparison of the results in 
the linear regime. Therefore, here the per particle cooling 
rate A is determined by one of the two cooling curves 
adopted in Paper I. However, because we study nonlinear 
effects in this paper, we do not represent the cooling curves 
as power laws in the neighborhood of the equilibrium tem- 
perature, and instead we implement the full nonlinear curve 
over a wide range of temperature. Both of the cooling 
curves used in the simulations are plotted in Figure 1. Each 
simulation uses either the cooling curve for interstellar gas 
calculated by Dalgarno & McCray (1972; hereafter DM) 
plotted in Figure la or the curve described by MacDonald 
& Bailey (1981; hereafter MB) for photoionized gas of 
reduced metallicity plotted in Figure lb. In both cases, the 
heating rate H is then determined by the requirement that 
the jet be in equilibrium initially, i.e., the right-hand side of 
equation (3) is zero. 

2.3. Initial and Boundary Conditions 

The simulations presented in this paper are all computed 
using Cartesian coordinates. We vary the size of the compu- 
tational domain in the axial direction depending on the 
Mach number of the jet and the perturbation frequency in 
order to ensure that many wavelengths of the most unstable 


mode can be captured. Typically, the grid is of size 
— 10R jt < Y < 10jRj t in the transverse direction and 0 < 
X < 100 — 400P jt in the axial direction, where R jt is the 
radius of the jet. In different simulations, we use from 200 to 
600 zones in the transverse direction, with from 10 to 40 
zones across a jet diameter. To allow the maximum 
resolution within the jet beam, nonuniform-sized zones are 
used beyond | Y | > 2R jt . In the axial direction, 400-800 
zones are used. Outflow boundary conditions are used 
along the outer Y boundaries, along the outer X boundary, 
and along the inner X boundary for | Y \ > 1 R jt . For | Y | < 
lR jt along the inner X boundary, inflow boundary condi- 
tions are used with the variables held fixed at the values 
appropriate to the initial equilibrium structure of the jet. 

We initialize the jet across the entire mesh for | Y \ < R jt , 
i.e., an equilibrium jet as opposed to a propagating jet. 
Initially, the jet is assumed to be perfectly collimated with 
density n jt , velocity v it , and temperature Tj t . The back- 
ground medium is stationary, with a uniform density n ex 
and with a pressure that matches that of the jet so that 
T ex = (n jt /n ex )Tj t . Motivated by recent observations of 
various protostellar jets (see, e.g., Ray 1996), we adopt as 
typical values for a protostellar jet a radius R it = 2.5 x 10 15 
cm; external Mach number M ex = v jt /a ex = 5 or 20, where 
a ex is the sound speed in the external medium; number 
density n Jt = 600 cm -3 ; and temperature Tj t = 1000 K. The 
ambient interstellar gas has a number density n ex = 60 
cm -3 , temperature T = 10,000 K, and sound speed a ex = 
1.18 x 10 6 cm s _1 . This implies we are studying overdense 
jets, with n it /n ex =10. There are two important timescales in 
the simulations. The first is the dynamical or sound crossing 
time t = Rjt/flj,, which for the parameter value listed above 
is equal to 213 yr. The second is the grid crossing time 




Log ( T ) 


Log ( T ) 


Fig. 1.— Per particle cooling rates used in the simulations for the (a) Dalgarno & McCray (1972; DM) cooling curve, and ( b ) MacDonald & Bailey (1981 ; 
MB) cooling curve. 
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t c ~ L grid /y jt , where L grid is the axial length of the computa- 
tional domain. This latter timescale varies with the Mach 
number of the jet and the size of the grid; each simulation is 
followed for several t c so that the jet reaches a quasi-steady 
structure in the nonlinear regime. 

As discussed in Paper I, in this study we focus on the 
development of K-H modes in jets that are overdense with 
respect to their surroundings. Because we must study 
pressure-matched jets in order to begin with an equilibrium 
state, this implies the jet beam is cooler than its surround- 
ings. Physically, our models may represent a dense jet beam 
embedded in a hotter and less dense cocoon that has been 
swept up by the bow shock associated with the leading 
surface. Dense molecular outflows that are often observed 
to be associated with protostellar jets (Reipurth & Cerni- 
charo 1995) may surround the cocoon. There is substantial 
observational evidence that protostellar jets are overdense 
with respect to their surroundings (see, e.g., Ray 1996). The 
fact that they remain well collimated for lengths that are 
much larger than their radii argues they must be in pressure 
equilibrium with their surroundings (or they would expand 
rapidly). Finally, the kinematic model of a dense jet embed- 
ded in a less dense cocoon that is in turn ensheathed by a 
very dense and cold molecular gas is consistent with 
detailed observation of a number of systems, including HH 
111 (Nagar, Vogel, & Stone 1996). 

In order to excite K-H modes, at time t = 0, a sinusoidal 
perturbation of the transverse velocity is turned on at the 
nozzle (x = 0). (As a test of our numerical methods, we have 
verified that the code will hold the initial equilibrium state 
indefinitely if no perturbation is applied to the jet.) The 
amplitude of the perturbation is set to u y = 0 jt M ex a ex , 
where 6 n is varied from 0.0025 to 0.01 for different driving 
frequencies. The velocity perturbation is of the form 

u y (x = 0) = 0 jt M ex a ex sin {cot) . (6) 

Because the wavelength of the most unstable modes varies 
with perturbation frequency, the axial length of the grid is 
chosen to ensure the grid will encompass several unstable 
wavelengths. 

3. RESULTS 

3.1. Evolution in the Linear Regime: The Asymmetric 
Surface Wave 

A linear analysis of the stability of the adiabatic slab jet 
reveals that the dominant K-H instable waves in a thermal 
pressure-confined jet are the asymmetric surface and sym- 
metric and asymmetric internal body waves (see, e.g., 
Hardee & Norman 1988). In general, the same is true for the 
cooling slab jet; however, the growth rates and wavelengths 
of unstable modes are strongly affected by the effect of 
cooling (Paper I), and at high frequency, even the symmetric 
surface wave can be important on the cooling jet. As an 
important consistency check between the analytic calcu- 
lations presented in Paper I and the numerical simulations 
presented here, we have compared the growth rates of K-H 
modes computed using both methods in the linear regime. 
We concentrate on the asymmetric surface wave, since in 
the frequency domain we are studying, the surface wave 
dominates the internal body waves. Moreover, we study the 
response of the jet to different driving frequencies, i.e., we fix 
co and study spatial instabilities. While the temporal growth 
rate can contribute to the total spatial growth rate (see, e.g., 


Zhao et al. 1992), in the linear regime the contribution is 
negligible. 

We have found a straightforward and accurate way to 
measure the linear growth rate of waves in our hydrody- 
namical simulations based on evolution of the “ color ” vari- 
able C that labels the jet and ambient gas. Perturbations 
induced at the jet nozzle result in a transverse velocity field 
that can be written as 

u y (x 9 t) = 0 jt M ex a ex exp {kjx} sin (cot) , (7) 

where kj is the imaginary part of the wavenumber. For a 
uniform ambient gas, the displacement amplitude of the 
surface wave is therefore given by (see, e.g., eq. [11] in 
Hardee, Norman, & Clarke 1994) 

A = A 0 exp {kj x} . (8) 

A schematic illustration of the growth of the surface wave is 
given in Figure 2. The locus of points marking the 
maximum transverse displacement of the jet at each axial 
location should therefore follow the exponential curve 
described by equation (8). We use the color variable to 
measure this exponential envelope in the following manner. 
The jet is evolved for many dynamical times, and the 
maximum transverse amplitude that the jet material 
reached throughout the evolution (found by plotting the 
envelope that encloses material with C = 1) at each axial 
position is measured. We find that a plot of this amplitude 
versus axial distance from the nozzle follows the exponen- 
tial envelope illustrated in Figure 2 extremely well. Figure 3 
demonstrates this procedure for a simulation of a Mach 5 
cooling jet computed using the MB cooling curve and per- 
turbed at a frequency of <x>R jt /a ex = 0.5. 

In order to compare the analytic and numerical growth 
rates at a variety of frequencies and jet parameters, we have 
performed numerical simulations at coR n /a ex = 0.2, 0.5, 2.0, 
and 5.0 for both Mach 5 and Mach 20 jets computed with 
an adiabatic equation of state, DM cooling, and MB 
cooling. The parameters of each of these test simulations are 
listed in Tables 1 and 2. In these tables, an upper limit to the 
perturbation wavelength is given by /l max = 2nu/coR }t . The 
true perturbation wavelength depends on the wave propa- 
gation speed and is different from the above estimate by the 
factor vju. The linear analysis indicates that to zeroth- 
order a reasonable estimate of the wave propagation speed 
is v w = [rj 1/2 /(l + r} ll2 J]u , where rj = p jt /p ex . A more accu- 
rate determination requires referral to the complete disper- 
sion relation solutions given in Paper I. We choose 0 jt to be 
small, but since the displacement amplitude of the surface 



Fig. 2. — Illustration of the growth of the fundamental surface wave. 
The solid sinusoidal line denotes the perturbation applied to the jet beam. 
The dotted line represents the effect of the exponential growth of the 
amplitude of the perturbation. The dashed line denotes the envelope of the 
maximum amplitude of the transverse position of the jet beam at each axial 
location. 
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Fig. 3.— Illustration of our procedure for measuring the growth rate of 
K-H modes from the numerical simulations. Top: The dotted line shows 
the isocontour of C = 1 at a late time in a simulation of a Mach 5, DM 
cooling jet to demonstrate the typical structure of the jet beam. The bold 
line shows the maximum displacement of the jet beam (defined by the 
C = 1 isocontour) throughout the entire simulation. The solid line is an 
exponential fit to the maximum displacement. Bottom: Fit of the amplitude 
of the jet displacement. 


wave is proportional to Q it u/co, at high frequencies, 0 jt is 
increased in order to resolve the transverse displacement. 

In Figures 4 and 5, we plot the growth rates of the surface 
waves measured from the numerical simulations (plotted as 
points) with the analytic results of Paper I (plotted as lines) 
for the Mach 5 and Mach 20 jets, respectively. The top 
panel in each figure compares the analytical and numerical 
results for an adiabatic jet, the middle panel for a non- 
adiabatic jet computed with the DM cooling curve, and the 
bottom panel for a nonadiabatic jet computed with the MB 
cooling curve. 

From Figure 4, it can be seen that the growth rate of the 
sinusoidal mode surface wave measured from the numerical 
simulations agrees very well with the analytic growth rate 
calculated in Paper I for low frequencies ( coRjJu < 0.4; 
a>RjJa ex < 2.0). For the adiabatic and MB cooling jet, the 
growth rate of the sinusoidal mode body waves (plotted as 
dashed lines) exceeds that for the surface wave (plotted as a 
solid line) at high frequencies. Our numerically measured 

TABLE 1 

Parameters of Simulations of M = 5 Jets 
to toRjJu KJR i« * b /A., N >‘ Q j» 


0.2 0.04 157.0 2.5 10 0.01 

0.5 0.1 63.0 3.2 20 0.02 

2.0 0.4 15.7 3.8 40 0.03 

5.0 1.0 6.3 4.8 40 0.075 


a Upper limit to perturbation wavelength. 
b Length of the computational domain. 
c Number of zones across the jet diameter. 

TABLE 2 

Parameters of Simulations of M = 20 Jets 


0) 0 )RjJu ^max/^jt 


0.2 0.01 630.0 1.9 20 0.001 

0.5 0.025 250.0 2.0 20 0.0025 

2.0 0.1 63.0 6.3 20 0.005 

5.0 0.25 25.0 8.0 40 0.01 




coRj t /u 


Fig. 4. — Comparison of the analytic and numerical growth rates in the 
linear regime for a Mach 5 jet. In each panel, the solid line is the growth 
rate of the fundamental surface wave as calculated by the linear analysis in 
Paper I, while the dashed lines are the growth rates of the fundamental and 
first two body modes. Growth rates measured from the numerical simula- 
tions are plotted as points with error bars. Top : Comparison for an adia- 
batic jet. Middle : Comparison for a cooling jet computed with the DM 
cooling curve. Bottom : Comparison for a cooling jet computed with the 
MB cooling curve. 


growth rate agrees well with the analytic growth rate of the 
body waves at high frequency (coRjJu = 1; coJR jt /a ex = 5), 
although there is considerable error in the determination of 
the growth rate. This error results from the growth of multi- 
ple body waves in the jet, which makes measurement of the 
growth rate of the fastest growing wave difficult using the 
techniques described above and in Figure 3. 

Inspection of Figure 5, which compares the numerical 
and analytical growth rates for waves in a Mach 20 jet, 
reveals similar trends. At low frequencies (wR s Ju < 0.1; 
coR)Ja ex < 2.0), the numerical results agree very well with 
the analytic growth rate for the sinusoidal mode surface 
wave. At a higher frequency, coR }t /u = 0.25 {a>R jt /a ex = 5.0), 
the growth rate of the sinusoidal mode body waves exceeds 
that of the surface wave for the adiabatic and MB cooling 
jet, and our numerically measured rates agree well with the 
body wave growth rate. Note that by comparing Figures 4 
and 5, it is evident that the growth rates, cutoff frequencies, 
and frequency of maximum growth all scale as the ratio of 
the Mach number between the Mach 5 and 20 jets as 
expected (Paper I). 

These results confirm the linear prediction that the jet 
cannot respond to a driving frequency above about 
a>RjJa ex « 2 by the development of a large-amplitude 
surface wave distorting the entire jet. Instead, internal 
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Fig. 5. — Same as Fig. 4 for a Mach 20 jet 


waves can be generated at the corresponding short wave- 
length (e.g., behavior of the asymmetric surface mode shown 
in Fig. 5 in Paper I). Thus, the jet is dominated by internal 
waves from the faster growing body modes (adiabatic and 
MB cooling jets) or by internal waves from the cooling 
enhanced faster growing surface wave (DM cooling jet). 

3.2. Evolution in the Nonlinear Regime 

The primary purpose of this paper is to use numerical 
methods to follow the growth of K-H modes in cooling jets 
into the nonlinear regime and to study the structure in the 
saturated state. Because of the exponential growth of waves, 
the linear regime will be short lived in real systems; thus, the 
structure and morphology of observed jets should be best 
fitted by nonlinear hydrodynamical calculations. In the fol- 
lowing subsections, we describe the properties of the nonlin- 
ear stage of the K-H instability observed in our simulations. 

3.2.1. The Asymmetric Surface Wave 

The linear analysis shows that the growth of the surface 
wave experiences a maximum at some intermediate fre- 
quency in both adiabatic and MB cooling jets (see Figs. 4 
and 5). In the nonlinear regime, mode-mode coupling can 
excite perturbations at this frequency even if the jet is driven 
at a different rate. Since this frequency has the highest 
growth rate, in the nonlinear regime, jets can be dominated 
by structures at the accompanying wavelength. This maxi- 
mally growing frequency is sometimes referred to as the 
“resonant” frequency (see, e.g., Ferrari, Turssoni, & Zanin- 
etti 1981 ; Hardee 1987). 

To study the nonlinear response of the surface wave and 
investigate perturbations induced by mode-mode coupling, 


we have performed low-resolution simulations (10 zones 
across a jet diameter) of Mach 5 adiabatic (Ad), MB cooling 
(MB), and DM cooling (DM) jets driven at a relatively high 
frequency of coR^/u = 1.0. In each simulation, the computa- 
tional domain has dimensions 0 < X < 40(XR jt , and 

— 20R jt < Y < 20R jt , and a numerical resolution of 
200 x 400 grid points. By using low numerical resolution, 
we artificially suppress the growth of body waves and there- 
fore enhance the role of the surface wave in determining jet 
structure. The result is shown in Figure 6 (Plate 4), which 
presents gray-scale images of the logarithm of the density in 
each jet at a dynamical time of £ = 100 t (note the grid cross- 
ing time for these simulations is t c = 80t). From Figure 4, it 
is clear that the adiabatic jet has a weak maximum in the 
growth rate at a>R/u « 0.3, the DM cooling jet lacks any 
resonant frequency, and the MB cooling jet has a strong 
maximum in the growth rate at a>R/u « 0.3. The structures 
observed in Figure 6 seem to reflect the relative strength of 
the maximum in the growth rate. For example, the adia- 
batic jet shows a weak sinusoidal pattern, the DM cooling 
jet has a complex and time-dependent structure that is not 
dominated by any one frequency, and the MB cooling jet 
develops a large amplitude sinusoidal oscillation that 
undergoes exponential growth from the nozzle until the jet 
is disrupted near r = 300 J R jt . The quasi-steady oscillation 
pattern has a wavelength X & 26R jt that corresponds to a 
linear perturbation frequency coR }t /u & 0.2 despite the fact 
that perturbation frequency is much higher. This result is 
consistent with growth near the resonant frequency 
coR^/u = 0.29 and wavelength X = 17.5K jt predicted by the 
linear analysis but suggests that longer wavelengths that 
grow at less than the maximum growth rate can come to 
dominate as they can attain larger amplitudes. These results 
once again confirm the linear prediction that the jet cannot 
respond to a driving frequency above about coR n /a ex ^ 2 by 
the development of a large-amplitude surface wave at the 
corresponding short wavelength. 

In order to investigate the internal dynamics of the jet 
beam and asymmetric surface wave in more detail, we have 
calculated the evolution of Mach 5 adiabatic, DM cooling, 
and MB cooling jets driven at a frequency coR jt /u = 0.4 
(near the maximum growth rate of the surface wave) and 
using a computational domain of size 0 < X < 60R n by 

— 5P jt < Y < 5P jt and 200 x 400 grid points (giving 40 
zones across a jet diameter). We plot gray-scale images of 
the logarithm of the density at dynamical time t = 15t in 
Figure 7 (Plate 5) (the grid crossing time is t c = 12t). In 
each case, internal structure associated with sinusoidal 
oscillation grows strongly, forming a sequence of opposing 
shocks in the jet. The perturbation frequency of coR H /u = 
0.4 lies below the resonant frequencies of the body waves 
and the linear analysis predicts asymmetric surface mode 
wavelengths given by A/R jt = 13.1 (Ad), 14.3 (DM), and 12.8 
(MB) compared to the observed 13.3 (Ad), 14.4 (DM), and 
14.8 (MB), and with the exception of the MB cooling jet, the 
agreement is excellent. The linear analysis predicts compa- 
rable growth rates for the adiabatic and MB cooling jets 
and a significantly larger growth rate for the DM cooling 
jet. Note that the DM cooling jet appears to disrupt at a 
lesser distance, ^45K jt , than the adiabatic and MB cooling 
jets. In the cooling jets, the postshock gas becomes very 
dense. The corrugations in the edge of the jet produce 
oblique shock spurs in the ambient gas oriented in a 
sequence of opposing steps. Interestingly, the knots them- 
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selves are subject to Rayleigh-Taylor instabilities and 
slowly break up at late times. We discuss knot formation in 
K-H unstable jets more fully in § 3.2.4. 

3.2.2. The Asymmetric Body Waves 

In order to study the nonlinear response of the body 
waves, we have performed high-resolution simulations (20 
zones across a jet diameter) of Mach 5 adiabatic, MB 
cooling, and DM cooling jets driven at a relatively high 
frequency of coR^/u = 1.0. In each simulation, the computa- 
tional domain has dimensions 0 < X < 200 R jt and 
— 5R jt < Y < 5 R jt and a numerical resolution of 100 x 800 
grid points. By using high numerical resolution, we can 
study the role of both body and surface waves in determin- 
ing jet structure. The result is shown in Figure 8 (Plate 6), 
which presents gray-scale images of the density in each jet at 
a dynamical time of t = 40t. Note the similarity with the 
gross behavior from the low-resolution simulations shown 
in Figure 6 (i.e., the nearly symmetric adiabatic jet, the 
complex structure of the DM cooling jet, and the 
developing sinusoidal oscillation of the MB cooling jet). 
The developing oscillation of the MB cooling jet has the 
wavelength, X « 26R jt , observed in the lower resolution 
simulation (Fig. 6) at much longer dynamical time. These 
results confirm both (1) the linear prediction that the jet 
cannot respond to a driving frequency above about 
(wk jt /a ex « 2 by the development of a large-amplitude 
surface wave at the corresponding short wavelength and (2) 
the validity of the low-resolution simulations. 

The linear analysis predicts that at high frequencies, the 
growth of body waves is very rapid (see Figs. 3 and 4). 
Indeed, these high-resolution simulations of jets perturbed 
at high frequency show that the internal waves become non- 
linear only a few wavelengths from the nozzle. At the 
driving frequency of ooR^/u = 1.0, the linear analysis pre- 
dicts that the wavelength of the asymmetric mode surface 
wave is X/R n = 5.62 (Ad), 4.95 (DM), and 5.56 (MB) com- 
pared to the observed 2/R jt « 5.9 (Ad), 6.3 (DM), and 6.1 
(MB) at < 50R jt . Note that the driving perturbation should 
excite the surface wave preferentially to the body waves, and 
this could overcome lower growth rates (cf. the adiabatic 
and MB cooling jet growth rates in Fig. 4); thus, the inter- 
nal waves at axial distances of < 50R jt could be due to the 
surface wave. Differentiating between short-wavelength 
surface and body waves requires study of the internal struc- 
ture. Beyond about 50R jt , the diagonal patterns in the 
images become more complex and indicate the presence of 
mixed surface and body waves. At these larger distances, we 
see evidence for ripples on the jet surface with separation 
(wavelength) in the range ~ 6.3R jt -4.7jR jt and internal 
oblique shock structures with separations (wavelengths) in 
the range ~ 3R jt -2R jt on the order of one-half to one-third 
the separation of the ripples. We note that from the linear 
analysis a potential wavelength range for the asymmetric 
mode first body wave is X/R it < 5.4 (Ad), 7.2 (DM), and 5.7 
(MB), where the upper limit corresponds to the wavelength 
at maximum growth. The driving frequency is at a wave- 
length about 15% shorter. The linear analysis would also 
predict a potential wavelength range for the asymmetric 
mode second body wave to be X/R jt > 3.17 (Ad), 3.55 (DM), 
and 3.26 (MB), where the lower limit corresponds to 
maximum growth. The driving frequency is at a wavelength 
about 35% longer. Interestingly, the observed separations 
are shorter than would be expected from the linear analysis 


although well within the shorter wavelength growth rate 
range of these body waves. The observed patterns suggest 
that wave-wave coupling between the surface and body 
waves is not capable of prediction via a linear analysis. 

At much lower frequencies, the body waves are predicted 
by the linear analysis to be stable (adiabatic jet) or damped 
(cooling jet). In Figure 9 (Plate 7), we plot gray-scale images 
of the density in Mach 20 adiabatic and DM cooling jets 
driven at a low frequency of coR^Ju = 0.025. The computa- 
tional domain in these simulations is of size 0 < X < 
500 R jt by — 10R jt < Y < 10R jt , and the numerical 
resolution is 200 x 400 grid points (giving 20 zones across a 
jet diameter). Each jet is shown at a time of t = 25r 
(equivalent to one grid crossing time). Note that to reveal 
the transverse structures in more detail, we have distorted 
the images by reducing the axial distance by a factor of 4 
relative to the transverse distance. When viewed with the 
correct aspect ratio, both jets appear essentially straight, 
although the oscillations achieve an amplitude on the order 
of the jet diameter. With the axial ratio as plotted, the 
images in Figure 9 should be representative of a Mach 5 jet, 
owing to Mach number scaling of growth rates and wave- 
lengths. 

Unlike jets perturbed at high frequency, jets perturbed at 
this low frequency have a long linear evolution phase. The 
long-wavelength oscillations have an observed X « 180R jt 
(Ad) and 200R jt (DM) — we have used the outermost oscil- 
lation to determine the wavelength — compared to the pre- 
dicted 180R jt (Ad) and 242R jt (DM). The DM cooling jet is 
predicted by the linear analysis to have a growth rate 
^30% larger than the adiabatic jet, and, in fact, the adia- 
batic jet remains relatively well collimated, whereas the DM 
cooling jet is disrupted and forms high-density knots associ- 
ated with working surfaces in the ambient gas. The dense 
knots result from strong cooling effects (see § 3.2.4). In addi- 
tion to these surface waves, both jets show evidence for 
internal oblique structures with wavelength in the range 
40R jt -50R jt at axial distances in the range 100R jt -200R jt . 
The linear analysis in Paper I shows that the first asym- 
metric body mode wave on the adiabatic jet would respond 
to this low-frequency perturbation with a wavelength of 
84R jt , considerably longer than observed in the simulation, 
and that the first asymmetric body wave on the DM cooling 
jet is damped at this low frequency. However, we note that 
the observed wavelengths fall within the wavelength range 
23R jt -84R jt on the adiabatic jet and 32R jt -63R jt on the 
cooling jet, where the short-wavelength limit is at the 
maximum growth rate and the long-wavelength limit is at 
the perturbation frequency or at the lowest frequency at 
which the wave is still unstable. Thus, we identify these 
internal structures with asymmetric body waves excited by 
the low-frequency perturbation. 

In Figure 10 (Plates 8-10), we plot gray-scale images of 
the logarithm of the density and the color variable for simu- 
lations of Mach 5 adiabatic, DM cooling, and MB cooling 
jets. In each case, the jets are perturbed at a relatively low 
frequency of ooR^/u = 0.1, and the simulations are extended 
until dynamical time t = 200t (or 2.5 crossing times) when 
the jets have all reached a quasi-stationary state. We use a 
high-resolution grid consisting of 600 x 800 zones (20 zones 
across jet diameter) to resolve internal jet structure, with a 
computational domain of size 0 < X < 400R jt by 
— 30 R jt < Y < 30R jt . The linear analysis in Paper I pre- 
dicts that this frequency is below the body wave resonant 
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frequencies on the adiabatic jet and that the body waves are 
damped at this frequency on the cooling jets. The jets all 
show a long linear evolution phase and an overall sinus- 
oidal pattern associated with the growth of surface waves. 
The long-wavelength oscillations have an observed X « 
50R jt (Ad), 52R jt (DM), and 51R jt (MB) compared to a cal- 
culated 46.7R jt (Ad), 59.3R jt (DM), and 49.1R jt . The surface 
wave grows to large amplitude at large distances from the 
nozzle and disrupts the jet. Shock spurs driven into the 
ambient gas by the sinusoidal jet oscillations are evident. 

In addition to the surface wave, there are oblique struc- 
tures within the jets. In the cooling jets, these structures are 
manifested primarily as dense knots or filaments oriented 
on either side of the jet in an alternating pattern (see, e.g., 
Fig. 10c at x ~ 100R jt ). Strong cooling effects are clearly 
observable in these simulations, especially for the jet with 
MB cooling, where the cooling rate for jet material is very 
high. The strong cooling in this case leads to condensation 
in the region in which the amplitude of the perturbation 
becomes significant. The most prominent of these structures 
that are at an axial distance less than or on the order of 
100R jt have an observed wavelength X « 17.3R jt (Ad), 
13.4R jt (DM), and 14.2R jt (MB). The linear analysis in 
Paper I shows that the first asymmetric body mode wave on 
the adiabatic jet would respond to this low-frequency per- 
turbation with a wavelength of 21R jt , not much longer than 
that observed in the simulation. However, the first asym- 
metric body wave on the DM and MB cooling jets is 
damped at this low frequency. However, we note that the 
observed wavelengths fall very close to the longest wave- 
length 15.3R jt on the DM cooling jet and 14.7R jt on the MB 
cooling jet at which the first asymmetric body wave is 
unstable. Thus, we identify these internal structures with 
asymmetric body waves excited by the low-frequency per- 
turbation. 

3 . 2 . 3 . Disruption of the Jet 

Once the amplitude of a surface wave at the maximally 
growing wavelength becomes on the order of the jet radius, 
the jet can no longer remain collimated and will disrupt. 
For wavelengths longer than the maximally growing wave- 
length, the amplitude can be correspondingly larger before 
the jet disrupts. This result is a consequence of the fact that 
disruption occurs when transverse motions in the reference 
frame of the wave approach the sound speed in the jet fluid. 
The distance a jet propagates before disruption depends on 
its Mach number (higher Mach number jets propagate 
farther), on the growth rate of the surface wave at the per- 
turbation frequency (higher growth rates produce quicker 
disruption), on the initial perturbation amplitude, and on 
the perturbation wavelength relative to the maximally 
growing wavelength. Disruption of the jet by the surface 
wave can be observed in each of the jets shown in Figures 6 
and 10 and also in the DM cooling jet in Figure 9. Beyond 
the point of disruption, the jet and ambient gas are strongly 
mixed, and shocked jet and ambient material are com- 
pressed into dense knots. 

The results of an extensive series of simulations at differ- 
ent perturbation frequencies and amplitudes indicate that 
jet disruption can be summarized as follows : 

1. For very low frequency perturbations, jet disruption 
occurs via the surface wave or low-order body modes, even 
though both have low growth rates. 


2. For intermediate-frequency perturbations, the jets are 
disrupted by the surface wave. 

3. For high-frequency perturbations, cooling jets are dis- 
rupted by strong body waves. On the other hand, body 
modes in the adiabatic jet tend to saturate at a finite ampli- 
tude and do not disrupt the jet 

3 . 2 . 4 . F ormation of Knots 

A common feature observed in our simulations of cooling 
jets is the formation of dense knots of cool gas. We find 
there are two separate mechanisms by which such knots can 
be formed. The first is the sweeping up of ambient gas by 
either (1) large-amplitude oscillations of the jet caused by 
growth of the surface wave or (2) disruption of the jet. This 
mechanism is similar in principle to the processes that form 
knots at the working surface of propagating protostellar 
jets (see, e.g., Blondin, Fryxell, & Konigl 1990; Stone & 
Norman 1994). The second mechanism is associated with 
the growth of body waves internal to the jet or in the linear 
growth phase of the sinusoidal surface wave that is accom- 
panied by waves internal to the jet. In the nonlinear regime, 
these waves form shocks in the jet, and strong cooling in the 
postshock gas can result in condensation of dense knots. 
Note that in this case, the knots are composed entirely of 
cooled jet material and are contained within the jet, whereas 
in the former process, the knots are composed of both jet 
and swept up ambient material. 

The Mach 5 cooling jets shown in Figure 10 provide 
evidence for the formation of knots via both the first and 
second mechanism. Growth of the sinusoidal oscillation 
results in large-amplitude displacement of the jet that then 
sweeps up ambient gas. When enough gas has accumulated, 
the jet starts to break up, and mixing of jet and ambient 
material follows. As a result of strong cooling, dense knots 
are formed in the mixing region (especially for MB cooling, 
where jet material has a much higher cooling rate). This 
affect is seen quite strongly in the Mach 20 DM cooling jet 
shown in Figure 9 as well. Both of the Mach 5 cooling jets 
and the Mach 20 DM cooling jet in these figures demon- 
strate knot formation via the second mechanisms as well. In 
these jets, we observe strong oblique internal waves associ- 
ated with the first asymmetric body wave. The strength of 
the internal wave grows rapidly, and dense knots are 
formed inside the jet (particularly in the MB cooling jet), 
but the jets still remain well collimated. The Mach 5 cooling 
jets shown in Figures 7 and 8 also provide evidence for knot 
formation via the second mechanism. These jets also 
contain strong oblique internal waves, and dense knots are 
formed inside the jet before the overall sinusoidal oscillation 
reaches sufficiently large amplitude to disrupt the colli- 
mated outflow. 

Note that the excitation state of gas in knots should be 
different depending on whether the knot material came 
from jet fluid or from a mixture of jet and ambient fluid. 
Owing to the very high growth rate of the surface wave for 
the DM cooling jet, we expect knots to contain a mixture of 
jet and ambient fluid as the jet disrupts very quickly. In 
most of the cases in which perturbations have intermediate 
frequencies, we expect that both of the mechanisms can be 
responsible for the formation of dense knots. 

3 . 2 . 5 . Mixing and Entrainment 

In the nonlinear regime, K-H instabilities can strongly 
alter the structure of a collimated jet and lead to mixing and 
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entrainment of ambient gas. Such mixing transfers momen- 
tum from the high-velocity jet material to the ambient 
medium, accelerating the latter and increasing the mass of 
outflowing material. In the case of protostellar jets, entrain- 
ment of high-velocity jet material into the ambient medium 
has been suggested as a mechanism for driving low-velocity 
molecular bipolar outflows, either by steady entrainment 
due to turbulence in a mixing layer between the jet and 
ambient gas (Stahler 1994) or by prompt entrainment in the 
“turbulent wakes” behind multiple bow shocks (Raga & 
Cabrit 1993). In our simulations, we can use the Lagrangian 
tracer or “ color ” variable C of the fluid to study quantitat- 
ively the mixing of the jet and ambient gases. 

In Figure 1 la, we plot the entrainment volume V € (defined 
as the volume of fluid that has C > e, where 0 < e < 1 (see, 
e.g., Loken et al. 1996), while in Figure 11 6, we plot the 
entrained mass M d (defined as the mass of fluid with v x > <3, 
where v x is the axial component of velocity (see, e.g., Loken 
et al. 1996; Bassett & Woodward 1995) as a function of 
position x along the jet. Both plots are constructed at late 
dynamical times (t = 200t, or 2.5 grid crossing times) for a 
Mach 5 jet computed with DM cooling and driven at a 
frequency of coR jt /u = 0.1 (this is the jet shown in Fig. 106). 
The jet structure has reached a quasi steady state by this 
time, and we average the axial profiles of both variables 
over a short time interval to remove the influence of discrete 
knots. The plots in Figure 11 illustrate the spatial variation 
in entrainment along the length of the jet driven by the 
growth of K-H modes. 

As expected, we find that little mixing occurs near the jet 
nozzle, i.e., inside x = 100R jt , where the waves are still in 
the linear regime. Beyond 100R Jt , Figure 11a demonstrates 
that the entrainment volume V e shows a rapid increase for 
every value of e, which inidcates that strong mixing is 
occurring. We find that mixing of the ambient gas starts 
where the perturbation reaches significant amplitude and 
then grows exponentially with the amplitude of the pertur- 
bation. This spatial growth of the mixing fraction appears 
to be in agreement with the results of Bodo et al. (1995). The 
spikes in the entrainment volume occurring at x « 160 and 
x« 180 are due to discrete knots of fluid shed by the jet 
beam during the initial growth of the instability; they are 


not part of the main beam at late times (see Fig. 106). The 
entrained mass M s plotted in Figure 116 drops for v x > 5 
(the initial jet velocity), while M s increases for all values of 
v x < 5, and this clearly indicates mass is being entrained 
and accelerated by the jet. For very small values of v x (e.g., 
v x < 0.1), we find M d begins to increase within 50R jt of the 
nozzle, even before the jet disrupts and mixes significantly 
(i.e., before V e begins to increase). However, Figure 106 
shows that in this region, strong shock spurs are driven into 
the ambient gas by the sinusoidal distortion of the jet. These 
shock spurs result in acceleration of the ambient gas. 
Although most previous discussions of the entrainment and 
acceleration of ambient gas along a jet beam have focused 
on turbulence either in a mixing layer, or in the turbulent 
“ wake ” of internal bow shocks in pulsed jets, our simula- 
tions show that shock trains in the ambient gas can result in 
acceleration of ambient gas. These shock trains may be 
produced by the K-H instability, as our simulations demon- 
strate, or they may be the wings of internal bow shocks 
produced by pulsing the jet. In fact, all of our simulations of 
the asymmetric mode of the K-H instability show the for- 
mation of shock spurs, which accelerates the ambient gas to 
some degree. Moreover, we note that neither our hydrody- 
namical simulations of these shock spurs nor simulations of 
the wings of internal bow shocks (Stone & Norman 1994; 
Biro & Raga 1994) show evidence for “ turbulent wakes.” 

On the other hand, we find that strong entrainment of 
ambient gas does not occur until the jet disrupts. Moreover, 
since the jets are not always disrupted by the asymmetric 
modes, not all of our simulations show strong mixing and 
entrainment of ambient gas, e.g., the Mach 5 cooling jets 
shown in Figure 7. However, our two-dimensional simula- 
tions cannot capture all modes of the K-H instability 
present in three dimensions; e.g., we lose the higher order 
fluting modes. It is possible, therefore, that the mixing rate 
in the nonlinear regime will be higher in three-dimensional 
simulations of unstable jets. In addition, very high 
resolution two-dimensional simulations of the temporal 
growth of K-H modes in low Mach number adiabatic jets 
(Bassett & Woodward 1995) indicate that high-order modes 
can significantly increase the mixing rate, which suggests 
that numerical resolution may limit our results. However, it 




Fig. 11a Fig. 11 b 

Fig. 11.— (a) Entrainment volume V e vs. axial position at late times in a Mach 5 DM cooling jet perturbed at a frequency of ojR jt /u = 0.1 (cf. Fig. 10b). 
Curves shown are for e = 0.01, 0.02, 0.05, 0.1, and 0.5. (b) Entrained mass M 6 vs. axial position for the same jet. Curves shown are for S = 0.5, 1, 2, and 5. 
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is not clear if these results apply to the spatial growth of 
modes in high Mach number jets as studied here. 

4. COMPARISON TO PROTOSTELLAR JETS 

Although we have studied the dynamical effect of a 
cooling curve appropriate to photoionized gas of low metal- 
licity (appropriate for jets from active galactic nuclei such as 
Seyferts) in addition to a cooling curve more appropriate to 
protostellar jets, the most detailed observations of optical 
jets collected to date have been of protostellar systems. 
Thus, in this paper, we focused on simulations that use 
parameter values (see § 2.3) representative of protostellar 
rather than Seyfert jets. For example, we have studied over- 
dense (n jt /n ex =10) jets, whereas jets from active galactic 
nuclei are usually thought to be density matched or under- 
dense with respect to their surroundings. While detailed 
modeling of Seyfert jets is a fruitful direction for future 
work, in this section we focus on a comparison between our 
numerical simulations and observations of protostellar jets. 

There are a number of striking similarities between the 
structures we observe in our simulations and the morphol- 
ogy of observed systems. For example, the low-amplitude 
oscillations, the formation of asymmetric knots in the jet, 
and the shock spurs driven into the ambient gas resulting 
from the growth of the body waves (see Fig. 7) are all remi- 
niscent of features seen in HST observations of jets such as 
HH 34 and HH 111 (Reipurth et al. 1997). HST obser- 
vations of HH 47 reported by Heathcote et al. (1996) indi- 
cate a complex structure, with the jet undergoing several 
bends before terminating at the Mach disk and bow shock 
structure HH 47 A. Interestingly, the jet structure seen in the 
DM cooling jet simulation shown in Figure 9 strongly 
resembles the observations of HH 47. Fabry-Perot obser- 
vations of the kinematics of HH 47 (Hartigan et al. 1993) 
show that the fastest moving jet material is confined to the 
core of the beam, with slower moving material forming a 
thick sheath around the jet. Raymond et al. (1994) have 
compared the kinematics of the HH 47 jet to models based 
on both viscous entrainment and acceleration due to mass 
ejected from the beam by internal bow shocks. In this paper, 
we have found that shock spurs produced by low-amplitude 
kinks in the beam that result from the K-H instability may 
be another mechanism for accelerating ambient gas in a 
sheath surrounding the jet beam. Our simulations indicate 
that it is not necessary for protostellar jets to be 
“precessing” through a large angle to generate kinks or 
wiggles. Instead, small-amplitude linear waves introduced 
into the jet via small perturbations can excite K-H insta- 
bilities to produce similar structures. 

There is strong evidence, based on kinematic data and 
the bow shock-like shape, that many internal knots in pro- 
tostellar jets are due to pulsing. At the same time, many 
sources show low-amplitude wiggles and asymmetric knots 
(as well as internal bow shocks) that are difficult to explain 
by pulsing alone. Such jets may best be modeled by the 
growth of instabilities in a pulsed jet, a problem we have not 
considered here. 

The peak in the growth rate for sinusoidal surface waves 
on adiabatic and jets with a MacDonald & Bailey cooling 
curve is at a resonant frequency of o>R jt /a ex » 1.5. For jets 
with a Dalgarno & McCray cooling curve, this frequency is 
also about the highest frequency to which the entire jet can 
respond dynamically. A similar result is found for helical 
twisting in three dimensions. For parameters typical of pro- 


tostellar jets, the frequency coR jt /a ex « 1.5 corresponds to 
perturbations with coRJu ^0.1 and a period of approx- 
imately 200 yr. However, if the jet is driven from the inner 
regions of the accretion disk (r < 1 AU), the characteristic 
time for perturbations (of order an orbital time) are much 
shorter. Our simulations show that jets that are driven at 
frequencies much higher than the resonant frequency 
cannot respond dynamically with a short wavelength corre- 
sponding to the high frequency, although the jet may show 
irregular oscillation (see the DM cooling jet in Fig. 6) and 
form internal structures via internal waves (see Fig. 7). This 
result is consistent with the observation that most proto- 
stellar jets are remarkably straight and usually show no 
evidence of large scale sinusoidal oscillations (which in three 
dimensions would result from helical twisting of the jet) 
leading to disruption of the type shown in Figure 6. 

Of course, observed jets may contain structures associ- 
ated with three-dimensional modes of the K-H instability 
that cannot be captured in our two-dimensional simula- 
tions. We will report on the comparison of three- 
dimensional simulations of K-H unstable protostellar jets 
with HST observations in a future communication. 

5. CONCLUSIONS 

In this, the second of two papers, we have studied the 
asymmetric modes of the K-H instability in jets in which 
optically thin radiative cooling is dynamically important, 
such as the jets associated with protostellar outflows and 
Seyfert galaxies. In Paper I, we presented solutions to the 
dispersion relation in the linear regime for K-H modes in a 
cooling jet over a wide range of perturbation frequencies. 
Our results indicated that the linear growth rates of waves 
in a cooling jet are substantially different from those of an 
adiabatic jet. Depending on the details of the cooling func- 
tion adopted, the growth rate of waves could be either 
enhanced or reduced relative to those on the adiabatic jet. 
In this paper, we have continued this analysis by using 
time-dependent hydrodynamical simulations to follow the 
growth of unstable waves into the nonlinear regime. 

As a consistency check, we first compared the linear 
growth rates of waves in our numerical simulations with 
predictions from the analytical results of Paper I. We found 
excellent agreement between the measured and predicted 
growth rates over a wide range of frequencies. This obser- 
vation serves both as a code test, and a confirmation of the 
results of Paper I. 

In the nonlinear regime, we observe a variety of inter- 
esting effects. Growth accompanying low-frequency and 
long-wavelength perturbations can dominate the structure 
of the jet at late times, producing large-amplitude sinusoidal 
oscillations of the jet beam. Once the amplitude of the oscil- 
lations exceeds several jet radii, the jet is disrupted. On the 
other hand, higher frequency perturbations can excite body 
waves that result in the formation of internal sinusoidal 
structures. Strong internal shocks produced by these body 
waves can form a pattern of dense knots that are staggered 
asymmetrically along the length of the jet. We find that 
internal shocks produced by relatively short wavelength 
body waves can coexist on a jet with a much longer wave- 
length, slowly growing, sinusoidal oscillation of the entire 
jet. We also find that small-amplitude oscillations of the jet 
associated with a growing sinusoidal surface or with a 
strong body wave can drive shock spurs into the ambient 
gas, which leads to an alternating pattern of knots and 
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shock spurs. These shock spurs can accelerate the ambient 
gas without mixing it with the jet fluid. On the other hand, 
disruption of the jet by the K-H instability produces strong 
mixing and entrainment of ambient gas and the formation 
of dense knots of shocked ambient gas swept up by jet 
material. 

Structures observed in our simulations such as low- 
amplitude wiggles in the jet, dense knots formed asym- 
metrically with respect to the center of the jet, and an 
alternating pattern of shock spurs driven into the ambient 
gas are all reminiscent of high-resolution observations of 
protostellar jets. On the other hand, for a typical protostel- 
lar jet, dynamical periods associated with the inner (r < 1 
AU) regions of an accretion disk are much less than the 
period of perturbations to which the entire jet can respond 
dynamically (which is about 200 yr). We therefore expect 


protostellar jets to be dominated by the growth of body 
waves and internal asymmetric structures rather than the 
surface waves. It is the surface waves that lead to large-scale 
apparent sinusoidal oscillation (resulting from helical 
twisting) of the entire jet and that would result in disruption 
of the jet at nonlinear amplitudes, features that typically are 
not observed in real systems. 

In this initial study, we have confined our simulations to 
the two-dimensional slab jet. We will report on the linear 
stability analysis and hydrodynamic simulations of three- 
dimensional cooling jets in a future communication. 
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Fig. 6. — Gray-scale images of the logarithm of the density for a Mach 5 adiabatic (top), DM cooling (middle), and MB cooling (bottom) jet perturbed with 
a frequency wR-Ju = 1.0 at the nozzle. Note the MB cooling jet shows a regular sinusoidal pattern indicating the dominance of a resonance frequency as 
expected from the dispersion relation of the surface wave. 
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F IG 7 -Gray-scale images of the logarithm of the density for a Mach 5 adiabatic (top), DM cooling (middle) and MB cooling (bottom ) i jet perturbed with 
a frequency mRJu = 0.4 at the nozzle. Dense knots and filaments are formed in a sinusoidal pattern owing to the nonlinear evolution of the surface wave. In 


the cooling jets ( bottom two panels ), these knots are R-T unstable and fragment at late times. 
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Fig. 8.— Gray-scale images of the density for high-resolution simulations of a Mach 5 adiabatic (top), DM cooling (middle), and MB cooling (bottom) jet 
perturbed with a frequency c oR } Ju = 1.0 at the nozzle, showing the growth and interaction of body modes and surface waves. 
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Fig. 9. — Gray-scale images of the density for a Mach 20 adiabatic (top) and DM cooling (bottom) jet perturbed with a frequency coR-Ju = 0.025 at the 
nozzle. Internal structure produced by the body modes is evident. 
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Fig. 10a 

Fig. 10. — Gray-scale images of the logarithm of the density and the color variable C for a Mach 5 (a) adiabatic, (b) DM cooling, and (e) MB cooling jet 
perturbed with a frequency coR it /u = 0.1 at the nozzle. Strong mixing is evident once the jet disrupts. 
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Fig. 10c 
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